Phase-amplitude coupling: processing and estimation

Standard Imports


In [1]:
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

import seaborn as sns
sns.set_style('white')

Load data


In [5]:
x = np.load('./data.npy')
Fs = 512.

Visualize data


In [7]:
plt.figure(figsize=(15,5))
plt.plot(np.arange(0,10+1/Fs,1/Fs),x,'k')
plt.xlabel('time (s)')


Out[7]:
<matplotlib.text.Text at 0x18040be0>

Broad bandpass raw signal (4-200Hz)


In [46]:
# Process signal
from scipy.signal import firwin, filtfilt
def mybandpass(x,cf,Fs, w = 3):
    numtaps = w * Fs / np.float(cf[0])
    if numtaps % 2 == 0:
        numtaps = numtaps + 1
    taps = firwin(numtaps, np.array(cf) / np.float(Fs) * 2, pass_zero=False)
    return filtfilt(taps,[1],x)

cf = (4,200)
ecogBhi = mybandpass(ecogB,cf,Fs)
ecogDhi = mybandpass(ecogD,cf,Fs)

Remove high frequency artifacts


In [52]:
from scipy.signal import butter
def _butterstimrmv(x, cf, bw, rate=1000, order=3):
    '''
    Notch Filter the time series x with a butterworth with center frequency cf
    and bandwidth bw
    '''
    nyq_rate = rate / 2.0
    f_range = [cf - bw / 2.0, cf + bw / 2.0]
    Wn = (f_range[0] / nyq_rate, f_range[1] / nyq_rate)
    b, a = butter(order, Wn, 'bandstop')
    return filtfilt(b, a, x)

In [213]:
stimfreq = [65,128,195]
stimbw = [2,30,1]
Nfilters = len(stimfreq)
for nf in range(Nfilters):
    if nf == 0:
        ecogB2 = _butterstimrmv(ecogBhi,stimfreq[nf],stimbw[nf],rate=Fs)
        ecogD2 = _butterstimrmv(ecogDhi,stimfreq[nf],stimbw[nf],rate=Fs)
    else:
        ecogB2 = _butterstimrmv(ecogB2,stimfreq[nf],stimbw[nf],rate=Fs)
        ecogD2 = _butterstimrmv(ecogD2,stimfreq[nf],stimbw[nf],rate=Fs)

Plot power spectra


In [214]:
from voytoys.scott.spec_tools import fftmed
Hzmed = 0
xlim = (0,205)
f, psdB = fftmed(ecogB, Fs=Fs, Hzmed = Hzmed)
f, psdD = fftmed(ecogD, Fs=Fs, Hzmed = Hzmed)

f, psdBhi = fftmed(ecogBhi, Fs=Fs, Hzmed = Hzmed)
f, psdDhi = fftmed(ecogDhi, Fs=Fs, Hzmed = Hzmed)

f, psdB2 = fftmed(ecogB2, Fs=Fs, Hzmed = Hzmed)
f, psdD2 = fftmed(ecogD2, Fs=Fs, Hzmed = Hzmed)

In [215]:
plt.figure(figsize=(15,10))
plt.subplot(3,1,1)
plt.plot(f,np.log10(psdB),'k',label='before')
plt.plot(f,np.log10(psdD),'r',label='DBS')
plt.legend(loc='best')
plt.ylabel('Raw')
plt.xlim(xlim)

plt.subplot(3,1,2)
plt.plot(f,np.log10(psdBhi),'k',label='before')
plt.plot(f,np.log10(psdDhi),'r',label='DBS')
plt.ylabel('Bandpass 4-200Hz')
plt.xlim(xlim)

plt.subplot(3,1,3)
plt.plot(f,np.log10(psdB2),'k',label='before')
plt.plot(f,np.log10(psdD2),'r',label='DBS')
plt.ylabel('High frequency\n artifacts removed')
plt.xlim(xlim)


Out[215]:
(0, 205)

In [229]:
Hzmed = 1
f, psdB2 = fftmed(ecogB2, Fs=Fs, Hzmed = Hzmed)
f, psdD2 = fftmed(ecogD2, Fs=Fs, Hzmed = Hzmed)

plt.figure(figsize=(15,6))
plt.plot(f,np.log10(psdB2),'k',label='before')
plt.plot(f,np.log10(psdD2),'r',label='DBS')
plt.ylabel('High frequency\n artifacts removed')
plt.xlim(xlim)


Out[229]:
(0, 205)

Visualize time series


In [216]:
t = np.arange(0,10+1/Fs,1/Fs)
trange = (0,1)
tidx = np.logical_and(t>=trange[0],t<trange[1])
plt.figure(figsize=(15,5))
plt.subplot(2,1,1)
plt.plot(t[tidx],ecogB[tidx]-np.mean(ecogB[tidx]),'k',label='raw')
plt.plot(t[tidx],ecogBhi[tidx],'r',label='4-200Hz')
plt.plot(t[tidx],ecogB2[tidx],'b',label='artifacts\nremoved')
plt.legend(loc='best')
plt.subplot(2,1,2)
plt.plot(t[tidx],ecogD[tidx]-np.mean(ecogD[tidx]),'k',label='raw')
plt.plot(t[tidx],ecogDhi[tidx],'r',label='4-200Hz')
plt.plot(t[tidx],ecogD2[tidx],'b',label='artifacts\nremoved')


Out[216]:
[<matplotlib.lines.Line2D at 0x539de518>]

PAC comodulograms


In [217]:
from pacpy.pac import comodulogram
from matplotlib import cm

# Parameters
fp = (6,40)
fa = (20,200)
dp = 4
da = 8
pac_method = 'mi_tort'

f_phases = np.arange(fp[0], fp[1], dp)
f_amps = np.arange(fa[0], fa[1], da)

In [218]:
# Calculate comods
comodB = comodulogram(ecogB, ecogB, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodD = comodulogram(ecogD, ecogD, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodBhi = comodulogram(ecogBhi, ecogBhi, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodDhi = comodulogram(ecogDhi, ecogDhi, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodB2 = comodulogram(ecogB2, ecogB2, fp, fa, dp, da,fs=Fs, pac_method=pac_method)
comodD2 = comodulogram(ecogD2, ecogD2, fp, fa, dp, da,fs=Fs, pac_method=pac_method)

In [219]:
climraw = (0,.01)
climhi = (0,.01)
clim2 = (0,.01)

plt.figure(figsize=(15,15))
plt.subplot(3, 3, 1)
plt.pcolor(f_phases, f_amps+da, comodB.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.title('before DBS', size=20)
plt.yticks(np.arange(50,250,50),size=20)
plt.xticks(np.arange(10,40,10),size=20)
plt.xlabel('Phase frequency (Hz)', size=20)
plt.ylabel('Amplitude frequency (Hz)', size=20)

plt.subplot(3, 3, 2)
plt.pcolor(f_phases, f_amps+da, comodD.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.title('DBS', size=20)

plt.subplot(3, 3, 3)
plt.pcolor(f_phases, f_amps+da, comodD.T - comodB.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)
plt.title('change', size=20)

plt.subplot(3, 3, 4)
plt.pcolor(f_phases, f_amps+da, comodBhi.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)

plt.subplot(3, 3, 5)
plt.pcolor(f_phases, f_amps+da, comodDhi.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)

plt.subplot(3, 3, 6)
plt.pcolor(f_phases, f_amps+da, comodDhi.T - comodBhi.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)

plt.subplot(3, 3, 7)
plt.pcolor(f_phases, f_amps+da, comodB2.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)

plt.subplot(3, 3, 8)
plt.pcolor(f_phases, f_amps+da, comodD2.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)

plt.subplot(3, 3, 9)
plt.pcolor(f_phases, f_amps+da, comodD2.T - comodB2.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
#plt.clim(climraw)

plt.tight_layout()


figname = 'Charlescomod'
plt.savefig('C:/Users/Scott/Google Drive/Voytek/Voytek Lab Shared resources/manuscripts/PD/feb2016fig/'+figname + '.png')


Visualize high gamma amplitude trace


In [220]:
from pacpy.pac import pa_series
f_lo = (13,30)
f_hi = (50,200)
phaB2, ampB2 = pa_series(ecogB2, ecogB2, f_lo, f_hi, fs=Fs)
phaD2, ampD2 = pa_series(ecogD2, ecogD2, f_lo, f_hi, fs=Fs)

def _edgeadd_paseries(amp, fosc, Fs, w = 3):
    """
    Undo the removal of edge artifacts done by pacpy in order to align
    the extrema with their amplitudes
    """
    Ntaps = np.floor(w * Fs / fosc[0])
    amp2 = np.zeros(len(amp)+2*Ntaps)
    amp2[Ntaps:-Ntaps] = amp
    return amp2

phaB2 = _edgeadd_paseries(phaB2,f_lo,Fs)
ampB2 = _edgeadd_paseries(ampB2,f_lo,Fs)
phaD2 = _edgeadd_paseries(phaD2,f_lo,Fs)
ampD2 = _edgeadd_paseries(ampD2,f_lo,Fs)


C:\Users\Scott\Anaconda2\lib\site-packages\ipykernel\__main__.py:13: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
C:\Users\Scott\Anaconda2\lib\site-packages\ipykernel\__main__.py:14: DeprecationWarning: using a non-integer number instead of an integer will result in an error in the future

In [221]:
trange = (0,5)
tidx = np.logical_and(t>=trange[0],t<trange[1])

plt.figure(figsize=(15,6))
plt.subplot(3,1,1)
plt.plot(t[tidx],ecogB2[tidx])
plt.ylabel('Raw voltage')
plt.subplot(3,1,2)
plt.plot(t[tidx],ampB2[tidx])
plt.ylabel('HFA amp')
plt.subplot(3,1,3)
plt.plot(t[tidx],phaB2[tidx])
plt.ylabel('Beta phase')


plt.figure(figsize=(15,6))
plt.subplot(3,1,1)
plt.plot(t[tidx],ecogD2[tidx])
plt.ylabel('Raw voltage')
plt.subplot(3,1,2)
plt.plot(t[tidx],ampD2[tidx])
plt.ylabel('HFA amp')
plt.subplot(3,1,3)
plt.plot(t[tidx],phaD2[tidx])
plt.ylabel('Beta phase')


Out[221]:
<matplotlib.text.Text at 0x52e7c3c8>

Test comod


In [245]:
# Simulate a low and high frequency oscillation.
def norm01(x):
    return (x - np.min(x))/(np.max(x)-np.min(x))
from voytoys.scott.tools_sim import simphase
T = 2
Fs = 1000.
t = np.arange(0,T,1/Fs)
betaphase, beta = simphase(T, f_lo, dt=1/Fs, returnwave=True)
_, gamma = simphase(T, f_hi, dt=1/Fs, returnwave=True)
pacsig = norm01(beta) + gamma*norm01(np.cos(betaphase))*.2

plt.figure(figsize=(15,5))
plt.plot(t,pacsig)
plt.xlim((0,1))


fpac,psdpac = fftmed(pacsig, Fs=Fs, Hzmed = 0)
plt.figure()
plt.plot(fpac,np.log10(psdpac))
plt.xlim((0,250))


Out[245]:
(0, 250)

In [246]:
comodpac = comodulogram(pacsig,pacsig, fp, fa, dp, da,fs=Fs, pac_method=pac_method)

In [248]:
plt.figure(figsize=(7,7))
plt.pcolor(f_phases, f_amps+da, comodpac.T, cmap=cm.jet)
plt.axis([f_phases[0], f_phases[-1], f_amps[0]+da, f_amps[-1]])
plt.colorbar()
plt.yticks(np.arange(50,250,50),size=20)
plt.xticks(np.arange(10,40,10),size=20)
plt.xlabel('Phase frequency (Hz)', size=20)
plt.ylabel('Amplitude frequency (Hz)', size=20)


Out[248]:
<matplotlib.text.Text at 0x5e17f9e8>

2. Calculate PAC and shape metrics


In [104]:
# Calculate PAC
from pacpy.pac import ozkurt
pacB = ozkurt(ecogB2, ecogB2, f_lo, f_hi, fs=Fs)
pacD = ozkurt(ecogD2, ecogD2, f_lo, f_hi, fs=Fs)

In [119]:
print pacB, pacD


0.0235779570366 0.0252327200944

In [111]:
# Identify peaks and troughs
from voytoys.shape.shape import findpt, Esharp, ESRsharp, RDsteep
boundaryS = 100

pkB, trB = findpt(ecogB2, f_lo, f_lopass = None, Fs = Fs, boundary = boundaryS)
pkD, trD = findpt(ecogD2, f_lo, f_lopass = None, Fs = Fs, boundary = boundaryS)

In [114]:
# Calculate sharpness of extrema
ampPC = 0
spPB = Esharp(ecogB2, pkB, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)
spTB = Esharp(ecogB2, trB, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)
spPD = Esharp(ecogD2, pkD, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)
spTD = Esharp(ecogD2, trD, width,ampPC=ampPC, Fs=Fs, fosc = f_lo)

In [109]:
# Calculate extrema sharpness ratio on multiple time scales
esrmethod = 'aggregate'
width = 5

esrB = np.log10(ESRsharp(ecogB2,pkB,trB,width,esrmethod=esrmethod))
esrD = np.log10(ESRsharp(ecogD2,pkD,trD,width,esrmethod=esrmethod))

In [110]:
# Calculate rise & decay steepness
def RDsteepratio(Rsteep,Dsteep,logrds = True):
    if logrds:
        return np.log10(np.max((np.mean(Rsteep)/np.mean(Dsteep),np.mean(Dsteep)/np.mean(Rsteep))))
    else:
        return np.max((np.mean(Rsteep)/np.mean(Dsteep),np.mean(Dsteep)/np.mean(Rsteep)))
    
RsteepB, DsteepB = RDsteep(ecogB2, pkB, trB)
RDsteeprB = RDsteepratio(RsteepB,DsteepB)

RsteepD, DsteepD = RDsteep(ecogD2, pkD, trD)
RDsteeprD = RDsteepratio(RsteepD,DsteepD)

Visualize shape


In [118]:
plt.figure(figsize=(5,8))
plt.subplot(2,1,1)
Nbins=np.arange(-2,2,.1)
plt.hist(np.log10(spPB),Nbins,color='k',label='peaks',alpha=0.5)
plt.hist(np.log10(spTB),Nbins,color='r',label='troughs',alpha=0.5)
plt.ylabel('Count',size=20)
plt.legend(loc='best')
ax1 = plt.gca()
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax1.get_yticklabels(), visible=False)

plt.subplot(2,1,2)
plt.hist(np.log10(spPD),Nbins,color='k',label='peaks',alpha=0.5)
plt.hist(np.log10(spTD),Nbins,color='r',label='troughs',alpha=0.5)
plt.xlabel('Log Extrema sharpness',size=20)
plt.xticks(np.arange(-1,2,1),size=20)
plt.ylabel('Count',size=20)
plt.legend(loc='best')
ax1 = plt.gca()
plt.setp(ax1.get_yticklabels(), visible=False)

print sp.stats.ttest_ind(spPB,spTB)
print sp.stats.ttest_ind(spPD,spTD)


Ttest_indResult(statistic=-3.0895924037915368, pvalue=0.002138186849668914)
Ttest_indResult(statistic=-2.1188220763461514, pvalue=0.034704787014510324)

In [ ]: